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Abstract 

The  inverse  scattering  problem  is  solved  within  the  distorted  wave  Born  approximation  (DWBA) 
using  a  generalized  time  reversal  based  formulation.  The  theory  is  applicable  to  arbitrary  (non-point) 
Na  element  transmitter  arrays  and  Np  element  receiver  arrays  and  to  scattering  systems  characterized  by 
scattering  potentials  interrogated  in  a  set  of  scattering  experiments  employing  single  or  multiple  temporal 
frequencies.  The  inversion  scheme  generates  a  pseudo-inverse  reconstruction  of  the  scattering  potential 
based  on  the  SVD  of  the  Np  x  Na  multistatic  data  matrix  considered  as  a  mapping  from  the  Hilbert  space 
of  scattering  potentials  to  the  vector  space  of  Np  x  Na  x  Nu  complex  N-tuples,  where  JVW  is  equal  to  the 
number  of  discrete  temporal  frequencies  employed  in  the  set  of  scattering  experiments.  The  theory  and 
results  in  this  report  have  been  previously  published  and  appear  in:  1.  A.J.  Devaney  and  M.  Dennison, 
“Inverse  Scattering  In  Inhomogeneous  Background  Media,”  Inverse  Problems  19,  No.  4  (August  2003), 
855-870,  2.  M.  Dennison  and  A.J.  Devaney,  ‘Inverse  scattering  in  inhomogeneous  background  media:Part 
II  SVD  analysis  and  multiple  frequency  case,  Inverse  Problems ,  20,  2004,  1307-1324. 


1  Introduction 

The  research  reported  herein  was  conducted  over  the  past  three  years  and  consisted  of  the  extension  of 
the  methods  of  time  reversal  based  imaging  to  applications  of  inverse  scattering.  The  foundations  of  the 
developed  theory  were  published  in  2003  in  [1]  and  applied  to  the  inverse  scattering  problem  for  frequency 
dependent,  compactly  supported  scattering  potentials  formulated  within  the  distorted  wave  Born  approxima¬ 
tion  (DWBA)  [2].  In  that  paper  the  scattering  object  was  assumed  to  be  embedded  in  a  known  background 
medium  and  the  distorted  wave  Born  approximation  to  the  Lippmann  Schwinger  equation  was  shown  to 
lead  to  a  simple  expression  for  the  so-called  multistatic  data  matrix  K (/3; ,  a*,,  u).  This  latter  quantity  is  the 
scattered  field  amplitude  at  receiver  point  f3j  due  to  point  source  excitation  at  transmitter  point  a*  and 
temporal  frequency  w  and  was  assumed  to  be  known  (measured)  over  some  arbitrary  set  of  transmitter  and 
receiver  points  indexed  by  k  =  1, 2,  -  •  •  ,  Na,  j  =  1, 2,  •  •  • , Np  and  at  some  fixed  frequency  uj.  A  Hilbert  space 
formulation  of  the  inverse  scattering  problem  thus  stated  was  then  employed  that  lead  in  a  simple  way  to 
an  inversion  algorithm  that  generates  a  least  squares,  pseudo-inverse  of  the  scattering  potential  at  the  given 
temporal  frequency  w. 

The  theory  and  results  established  in  [1]  were  extended  in  [3]  to  be  applicable  to  non-point  transmitter 
and  receiver  elements  and  to  dispersionless  scattering  potentials  for  which  the  scattering  data  are  known 
at  multiple  temporal  frequencies  uif,  f  =  1,2,--  -  ,NU.  Although  both  generalizations  were  found  to  be 
simple  extensions  of  the  underlying  theory  presented  in  the  earlier  paper  and  require  only  minor  changes  in 
notation  and  implementation  they  can  be  important  in  practical  applications  that  generally  employ  extended 
2D  antenna  or  transducer  elements  and  broadband  or  stepped  frequency  excitation. 

The  least  squares,  pseudo-inverse  inversion  schemes  developed  in  the  grant  are  based  on  the  singular  value 
decomposition  (SVD)  of  the  underlying  DWBA  mapping  V  :  Hy  — >  CN,  where  Hy  is  the  Hilbert  space  of 
compactly  supported  scattering  potentials  and  CN  the  space  of  complex  N  tuples  where  N  =  Np  x  Na  x  Nu. 
I  emphasize  that  the  use  of  the  SVD  employed  here  to  solve  the  inverse  scattering  problem  is  vastly  different 
from  its  conventional  use  in  other  time-reversal  based  imaging  and  inverse  scattering  applications  [4,  5,  6] 
where  it  is  used  to  decompose  the  linear  mapping  from  the  incident  field  to  the  scattered  field.  In  the  study 
reported  here,  it  is  used  to  decompose  the  mapping  from  the  Hilbert  space  of  compactly  supported  scattering 
potentials  to  the  finite-dimensional  vector  space  representing  the  scattered  field  data.  These  two  mappings, 
and  associated  SVD’s,  are  totally  different  and  are  directly  applicable  to  different  classes  of  problems. 
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2  Problem  Formulation  and  Review  of  [1] 

I  consider  the  general  scenario  of  a  scattering  object  compactly  supported  in  some  volume  V  that  is  embedded 
in  a  known  inhomogeneous  background  medium  and  is  interrogated  in  a  set  of  Na  scattering  experiments 
where  in  each  experiment  the  incident  wave  is  generated  by  a  point  source  located  at  ak,  k  =  1, 2,  •  •  •  ,  Na 
and  the  scattered  wave  is  measured  over  a  set  of  Np  point  receiving  sensors  located  at  /3j,  j  =  1, 2,  •  •  • ,  Np. 
The  suite  of  data  obtained  from  these  experiments  is  called  the  multistatic  data  matrix  and,  within  the 
so-called  distorted  wave  Born  approximation  (DWBA)[2],  is  related  to  the  object’s  scattering  potential  via 
the  equation  [1] 

Ay.fcM  =  [  d3rG0(/3j,r,u!)O(r,u)G0(r,ak,u}),  (1) 

Jv 

where  u>  is  the  temporal  frequency,  0(r,  u)  is  the  object’s  (generally  complex)  scattering  potential  relative  to 
the  known  background  and  Gq  is  the  background  Green  function  which  is  assumed  known.  The  multistatic 
data  matrix  Kj<k  can  be  directly  measured  at  one  or  more  frequencies  to  via  monochromatic  source  excitation 
or,  alternatively,  can  be  computed  over  a  large  set  of  frequencies  via  Fourier  transformation  of  the  measured 
responses  from  pulse  excitation. 

By  defining 


Dn(w)  =  Kitk M,  n  =  j  +  (k-l)Np,  j  =  1, 2,  •  •  •  ,  Np,  k  =  1, 2,  •  •  • ,  Na,  (2) 

we  can  regard  the  relationship  Eq.(l)  to  be  a  linear  mapping  from  the  Hilbert  space  Hv  of  complex  valued 
scattering  potentials  0(r,  w)  compactly  supported  in  V  to  the  finite  dimensional  vector  space  CN(ui)  of 
complex  N-tuples  D  =  {Dn(iu)},  n  =  1, 2,  •  •  • ,  TV  where  N  =  Np  x  Na.  Following  the  treatment  presented 


in  [1]  I  express 

where  V  is  the 

Eq.(l)  in  the  form 

VO{oj)  =  D{u) 

integral  operator  Hy  — >  CN(w): 

(3) 

IT 

■b- 

*  e 

s- 

% 

II 

(4a) 

with 

/  \  j Go{Pj,r,aj)GQ(ak,r,u})  if  r  e  V 

7rn(r,a»)  —  \  , 

10  else- wise. 

(4b) 

The  above  set  of  equations  provide  the  underlying  formulation  of  the  inversion  scheme  developed  in  the  grant 
and  reported  in  [1]  and  in  [3]. 


2.1  Pseudo-inverse  of  the  DWBA  mapping 

At  each  frequency  oj  the  distorted  wave  Born  approximation  (DWBA)  mapping  Eq.(3)  defines  a  projection 
of  the  unknown  scattering  potential  0(r,u>)  onto  the  subspace  Hy  =  Span{^r„)  C  Hxi  of  Hy  that  is  spanned 
by  the  set  of  functions  {7rn(r,  w)},  n  =  1,2,  -  -  -  ,TV.  This  fact  was  the  basis  for  the  pseudo-inverse  of  this 
mapping  presented  in  [1]  where  the  pseudo-inverse  0(r,ui)  was  expanded  into  the  set  of  spanning  functions 
7 rn  with  coefficients  Cn  that  were  computed  by  direct  matrix  inversion.  Thus,  in  particular,  we  have  that 

N 

d(r,w)  =  ^Cn(w)7r„(r,w)  (5a) 

n=l 


where 


and  where 


N 

<  7rn,0  > H\>—  Dn(cu)  =  )  ’  <  7Tn,  7T„'  Gn' 

n'= 1 


<  fi,f 2  >HV=  f  d3r  /i  (r)/2(r) 
Jv 


(5b) 

(6) 
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is  the  standard  inner  product  in  Hy.  Eq.(5b)  can  be  expressed  in  the  matrix  form 

n(w)C'(w)  =  D(u) 


(7) 


where  II (w)  =<  7rn,  7rn<  >hv  is  the  N  x  N  inner  product  matrix  of  the  spanning  functions  7rn  and  C(oj)  = 
[Ci (u),  C^oj),  -  -  -  ,  Cjv(w)]t  is  the  column  vector  formed  from  the  expansion  coefficients  of  the  pseudo-inverse 
Eq.(5a). 

The  final  step  of  the  inversion  consists  of  diagonalizing  the  inner  product  matrix  II  (w)  and  generating  a 
least  squares  inversion  of  Eq.(7);  viz., 

n  =  UAU*  =>  6  =  UA+U'D,  (8) 

where  A+  denotes  the  diagonal  matrix  whose  non-zero  diagonal  elements  are  the  reciprocals  of  the  non-zero 
eigenvalues  A„  >  0  of  the  II  matrix.  As  discussed  in  [1]  the  expansion  Eq.(5a)  with  expansion  coefficients 
determined  from  Eq.(8)  is  a  least  squares,  pseudo-inverse  of  the  inverse  scattering  problem  within  the  DWBA. 
That  O  is  a  pseudo-inverse  (minimum  L2  norm  solution)  follows  from  the  fact  that  the  subspace  Hu  is  the 
pre  image  space  (orthogonal  complement  of  the  null  space)  of  the  DWBA  mapping  V.  The  fact  that  it  is 
a  least  squares  solution  follows  from  the  use  of  a  least  squares  solution  to  the  matrix  equation  defining  the 
expansion  coefficients. 


3  Generalizations  and  Extensions 


The  inversion  scheme  developed  in  [1]  and  summarized  above  admits  a  number  of  generalizations  and  ex¬ 
tensions  as  reported  in  [3].  Here  we  will  only  discuss  the  extension  to  non-point  transmitter  and  receiver 
elements  and  to  the  multiple  frequency  case  for  dispersionless  scattering  potentials. 

The  basic  scheme  is  easily  generalized  to  be  applicable  to  experiments  employing  extended  planar 
(non-point)  transmitter  and  receiver  elements  upon  simple  replacement  of  the  background  Green  functions 
G0(/3j,r,u})  and  Go{r,ak,uj)  by  the  wave  fields 

So(Pj,r,u)-  J  d2r'1Zr{f3j,T' ,  u;)Go(r',  r,  w),  (9a) 

g0(r,ak,u>)  =  J  d2r'TZt{r',ak,u))Go{r,r\tj).  (9b) 

The  wave  field  Qo(/3j,r,uj)  corresponds  to  the  measured  response  at  an  extended  planar  sensor  element  lo¬ 
cated  at  (3j  and  having  receiving  response  function  TZr(/3j,  r',  lj)  while  Go{r,  Qt,w)  is  the  wave  field  generated 
by  an  extended  planar  source  located  at  ak  and  having  transmission  response  function  77((r', ak,uj).  For 
this  more  general  situation  the  multistatic  data  matrix  as  defined  in  Eq.(l)  has  to  be  replaced  by 

=  J  d2r"  J  dV  UriP^r",  ui)KT»:T'  (w)77t (r',  ak,  ui) 

=  [  d3rg0{/3j,r,uj)0{r,(j)go(r,ak,u ’) 

Jv 

where 

KT"y{bj)  —  j  d3r  Gq(t",  r,  w)C>(r,  u>)G0(r,  r',  w), 

Jv 

is  the  multistatic  data  matrix  evaluated  between  the  source  point  r'  and  receiver  point  r".  With  the 
replacement  of  K]  k  with  fCjik  and  the  spanning  functions  irn  by 


Mn(r,w) 


f55(/3j,r,w)g5(r,afc,w) 

1° 


if  r  €  V 
else-wise. 


the  entire  theory  is  applicable  to  experiments  employing  non-point  transmitter  and  receiver  elements  char¬ 
acterized  by  the  transmission  and  receiver  response  functions  TZt  and  TZr. 
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4  Multi-frequency  case  for  frequency  independent  potentials 

In  the  theory  developed  in  [1]  and  outlined  above  the  object  scattering  potential  is  assumed  to  be  dependent 
(in  an  unknown  way)  on  the  temporal  frequency  w  so  that  Eq.(3)  is  a  mapping  from  the  Hilbert  space  Hy 
of  compactly  supported  scattering  potentials  at  fixed  frequency  to  to  the  N-dimensional  vector  space  Cn(uj) 
of  complex  N-tuples  D(lj)  =  {Z)n(o;)}  at  that  specific  fixed  frequency.  In  the  special  case  where  the  object 
scattering  potential  is  independent  of  cj  then  the  basic  structure  of  this  mapping  changes  and  becomes  a 
mapping  from  the  Hilbert  space  Hy  of  frequency  independent  scattering  potentials  O(r)  compactly  supported 
in  V  to  the  tensor  product  of  the  N—  tuple  spaces  CN (uS)  over  the  band  of  frequencies  used  in  the  set  of 
experiments. 

We  can  restructure  the  basic  DWBA  mapping  within  this  new  formulation  as  follows.  First,  lets  assume 
that  we  have  experimental  data  at  discrete  frequencies  u>  =  w/,  /  =  1,2,  -,  Nw.  In  the  spirit  of 

our  definition  Eq.(2)  we  define 

Dn  =  Kjik(uf),  n  =  j  +  (k  -  l)Np  +  (f  -  1)N„,  (10) 

where,  as  before,  j  =  1, 2, •••  ,Np,  k  =  1,2, •••  ,Na.  Eq.(10)  defines  an  AT— tuple  with  AT  now  equal  to 
N  =  Np  x  Na  x  Nu  and  the  DWBA  mapping  for  frequency  independent  scattering  potentials  becomes 


VO  =  D 

(11) 

where  V  is  the  integral  operator 

V=  f  d3r<( r), 

Jv 

(12a) 

with  Trn(r)  defined  as 

vrn(r)  = 

1  r,  u >f)G’0{ak,  r,  uif) 

1° 

if  r  G  V 

else- wise. 

(12b) 

With  the  above  definitions  the  pseudo-inverse  employed  for  frequency  dependent  potentials  can  be  applied 
with  minor  modification.  Thus,  in  particular,  we  have  as  before 


N 

0( r)  =  ^2  cnitn{r)  (13a) 

71=1 

where  the  expansion  coefficients  Cn  are  solutions  to  the  coupled  set 

N 

Ai  =  )  [  Cn>  <  7Tn,7rn'  •  (13b) 

7l'  =  l 

Eq.  13b  can  be  re-written  in  matrix  form  as 

D  =  UC  (13c) 

where  n(n,  n')  —<  7rn,  irni  >hv  is  the  N  x  N  inner  product  matrix.  The  least  squares  pseudo-inverse 
solution  for  the  scattering  potential  is  then  given  by  Eq.(13a)  with  the  expansion  coefficients  given  by  the 
least  squares  solution  of  Eq.(13c). 

It  is  clear  from  the  above  that  from  a  strictly  computational  viewpoint  the  inversion  of  frequency  inde¬ 
pendent  scattering  potentials  is  structurally  and  formally  identical  to  that  of  frequency  dependent  potentials. 
The  effect  of  adding  scattering  data  acquired  at  different  frequencies  is  treated  in  the  same  way  as  would  be 
the  addition  of  scattering  data  acquired  at  additional  transmitter  or  receiver  points.  Both  of  these  result  in 
an  increase  in  the  number  of  expansion  functions  {7rn}  generating  the  subspace  Hu  =  Span{7r„}  of  Hy  and, 
presumably,  increase  the  dimension  of  this  subspace  and  the  quality  of  the  inversion.  Although  transmitter 
and  receiver  location  diversity  and  frequency  diversity  are  formally  equivalent  for  non-dispersive  scatter¬ 
ing  potentials,  they  are  not  equivalent  in  terms  of  information  content.  In  particular,  frequency  diversity 
is  generally  of  more  importance  in  reflection  type  experiments  such  as  are  used  in  pulse-echo  ultrasonics, 
GPR  (ground  penetrating  radar)  and  seismic  exploration  while  source  and  receiver  diversity  are  critical  in 
transmission  type  experiments  employed  in  tomographic  type  experiments  [8] . 
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5  Reformulation 


The  computer  implementation  of  the  single  or  multiple  frequency  inversion  algorithms  involves  the  following 
steps: 

1.  Computation  of  the  inner  product  matrix 

n(n,n')=  /  d3r<(r)jrn-(r) 

Jv 

where  n,n'  =  1, 2,  •  •  •  ,  N  with  N  =  Np  x  Na  x  Nu, 

2.  Diagonalization  of  the  inner  product  matrix:  II  =  UAU', 

3.  Computation  of  the  least  square,  pseudo-inverse  C  of  the  expansion  coefficient  column  vector  C:  C  = 
UA+WD, 

4.  Computation  of  the  least  squares  pseudo-inverse  O(r): 

N 

6(r)  =  CnVn(r). 
n— 1 

Steps  1  and  2  are  data  independent  and,  hence,  can  be  performed  prior  to  actual  data  acquisition  and, 
thus,  do  not  affect  the  execution  time  of  the  inversion  algorithm.  Steps  3  and  4,  on  the  other-hand,  are 
performed  after  or  in  conjunction  with  data  acquisition  and  can  impose  rather  severe  demands  both  in  terms 
of  execution  time  and  computer  storage  requirements.  Step  4,  in  particular,  requires  the  storage  of  the  set 
of  spanning  functions  irn{r),  n  =  1, 2,  •  •  • ,  N  and  Vr  G  V  which  can  require  a  great  deal  of  storage  space  for 
three-dimensional  (3D)  problems  and/or  in  the  multi-frequency  case  where  the  number  N  =  Np  x  Na  x  Nu 
can  be  quite  large  even  for  moderately  large  numbers  of  transmitter  and  receiving  elements.  These  high 
computational  costs  stem  from  the  fact  that  the  set  of  spanning  functions  {itn}  are  neither  orthogonal  or 
linearly  independent  with  the  result  that  the  inversion  is  not  economical  in  terms  of  computer  storage  or 
CPU  time.  This  problem  can  be  overcome  by  generating  a  linearly  independent  set  either  by  means  of  the 
Gram-  Schmidt  procedure  or  by  application  of  the  Singular  Value  Decomposition  (SVD)  to  the  basic  DWBA 
mapping  Eq.(ll).  The  SVD  approach  is  particularly  appropriate  and  will  be  employed  here  since  it  allows 
us  to  order  the  singular  functions  in  order  of  importance  according  to  the  size  of  their  associated  singular 
value.  Here,  we  will  employ  the  SVD  and  derive  a  new  formulation  of  the  pseudo-inverse  that  minimizes 
the  computational  burden  while  retaining  the  same  quality  of  inversion  as  obtained  by  the  basic  formulation 
presented  in  [l]and  outlined  above. 

5.1  SVD  of  the  DWBA  mapping 

To  compute  the  SVD  of  the  DWBA  mapping  it  is  necessary  to  introduce  an  inner  product  in  the  data  space 
CN  and  to  compute  the  adjoint  relative  to  the  two  inner  products  in  the  spaces  Hy  and  CN .  As  before 
we  employ  the  standard  inner  product  in  Hy  and  will  also  employ  a  standard  inner  product  in  CN;  i.e., 

<fuf2  >tfv=  f  d3r  ft  (r)/2(r),  V/j,  f2  €  Hv,  (14a) 

Jv 

N 

<  ui,v,2  >cN=  ^wI(n)u2(n),  Vui,u,2  €  CN.  (14b) 

n=l 

With  these  inner  products  we  find  that 

N  r 

<  U,Vf  >C*=  5ZU‘(n)  /  rf3r7rn(r)/(r) 

n=l  Jv 

r  N 

=  /  d3r  u(n)7r„(r)}*/(r)  =<Ptu,/>Hv 

JV  n=l 
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from  which  we  conclude  that 

N 

V^U  =  ^2  u(n)7Tn(r),  Vu  £  CN .  (15) 

71—1 

The  SVD  of  the  DWBA  mapping  employs  the  singular  system  { up  £  CN,vp  €  Hv,op  >  0},p  = 
1,2,--- ,  TV  where 

'Tv  p  —  cTpUp ,  V  ^  =  cTpVp.  (16a) 

Using  the  singular  system  we  then  find  that 

N 

V=Y1  aPuPvl’  (l6b) 

p=i 

where  vj,  :  H\>  — >  C  is  defined  in  the  usual  way  via  the  equation 

v\f  =  [  d3rVp(r)f(r)  =<  vp,f>Hv,  V/  £  Hv, 

Jv 

and  for  each  p  =  1,2,  ■  ■  ■  ,  N .  The  singular  vectors  vp  and  up  are  found  as  solution  to  the  normal  equations 

V^Vvp  =  ffpVp,  VV^up  =  ffpUp.  (16c) 

The  functions  {up(r),  av  >  0}  (with  appropriate  normalization)  form  an  orthonormal  basis  in  the  Hilbert 
space  Hv  while  the  N-tuples  (column  vectors)  {up,  ap  >  0}  form  an  orthonormal  basis  in  CN .  Moreover,  the 
set  {vp,  <jp  >  0}  provides  an  orthonormal  basis  in  Hn  C  Hy  and  the  set  { up ,  ap  >  0}  provides  an  orthonormal 
basis  in  the  image  of  Hn  under  V.  The  set  {vp,  av  >  0}  are,  thus,  precisely  the  set  of  orthonormal  basis 
functions  needed  to  replace  the  set  of  spanning  function  {irn,  n  =  1,2,  ■  ■  ■ ,  N}  now  employed  in  the  pseudo¬ 
inverse  solution  developed  in  [ljand  outlined  above. 

5.2  Computing  the  Singular  Vectors 

We  can  write  the  l.h.s.  of  the  normal  equation  for  the  singular  vectors  up  in  Eq.(16c)  in  the  form 

r  N 

VV^Up  =  /  d3r  7r*(r)  up{n') 7r„-(r) 

Jv  n'=l 

N  r 

=  53  W  d3r>  KW)^n’(r')}up(n)  =  Tlup 
n'=l  Jv 

where  up{n)  is  the  nth  component  of  the  up  singular  vector  and  II  =<  nn ,  nn>  >hv  is  the  TV  x  TV  inner 

product  matrix  of  the  spanning  functions  7T„.  We  thus  see  that  the  normal  equation  for  up  is  simply  the 

matrix  equation 

II  Up  =  CTpUp  (17) 

which  is  the  diagonalizing  equation  for  the  n  matrix.  In  other  words  the  singular  vectors  up  are  the  eigen¬ 
vectors  of  the  II  matrix  and  the  square  of  the  singular  values  ap  are  the  eigenvalues. 

Once  the  singular  vectors  up  are  computed  the  associated  singular  vectors  vp  having  non-zero  singular 
values  Op  >  0  are  computed  using  Eq.(16a);  i.e., 

N  ,  , 

(TpVp  =  PfUp  =>vp=J2  (18) 

n=l  °P 

where  we  have  used  the  definition  of  the  adjoint  operator  P*  given  in  Eq.(15).  It  is  important  to  note  that 
the  singular  vectors  up,vp,  op  >  0  are  data  independent  and  can  be  computed  prior  to  data  acquisition 
and,  hence,  their  computation  will  not  affect  the  execution  time  or  storage  requirements  of  the  reformulated 
inversion  algorithm. 
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5.3  Generating  the  Inversion 

As  discussed  previously,  the  main  advantage  of  the  new  formulation  is  that  the  set  {vp},  op  >  0  are  or¬ 
thonormal  and,  hence,  form  a  minimal  set  of  functions  into  which  the  pseudo-inverse  O  of  the  scattering 
potential  can  be  expanded.  Thus,  in  particular,  it  is  only  necessary  to  store  the  singular  vectors  vp(r)  for 
<tp  >  0  rather  than  the  entire  spanning  set  7rn(r),  n  =  1, 2,  •  •  •  ,  N.  Moreover,  as  we  will  show  below,  only 
projections  <  vp,  O  > nv  of  the  scattering  potential  associated  with  singular  vectors  whose  singular  values 
cTp  exceed  some  threshold  value  <xp  >  e  >  0  can  be  reliably  computed  in  the  presence  of  additive  noise  or 
measurement  error  so  that  only  those  singular  vectors  whose  singular  values  satisfy  this  condition  need  be 
retained  and  used  in  the  inversion.  The  threshold  value  £  is  a  regularization  parameter  that  can  be  selected 
to  balance  noise  tolerance  versus  accuracy  of  the  inversion  algorithm.  We  mention  that  the  new  formulation 
still  requires  that  the  II  matrix  be  computed  and  then  diagonalized  as  in  the  old  scheme.  After  that  the 
eigenvectors  and  eigenvalues  of  II  are  used  to  compute  the  singular  vectors  vp  using  Eq.(18).  The  inversion 
can  then  be  generated  directly  from  the  stored  singular  system  {vp,  up,  ap  >  e}  as  we  will  now  describe. 

To  compute  the  pseudo-inverse  we  first  represent  the  data  vector  D  in  the  basis  up;  i.e., 

N 

D  =  ^  <  up>  D  >qn  up.  (19a) 

p=i 

In  a  similar  manner  we  can  expand  the  unknown  object  profile  O  into  the  set  of  singular  vectors  vp: 


N 

0=^2<vp,0  >HV  Vp.  (19b) 

p=  i 

If  we  now  substitute  this  set  of  equations  into  the  DWBA  mapping  VO  —  D  we  find  that 

N  N 

^2ap<  Vp ,  O  >HV  Up  =  ^2<Up,D  >c jv  Up  (20a) 

p=i  p=i 


from  which  we  conclude  that 


<  Vp,  O  > fjv 


<  Up,  D  >cn 

(Tp 


if  op  >  0.  The  pseudo- inverse  is  thus  given  by 


°=E 

<Tp>0 


<  Up,  D  >cn 

Op 


vv. 


(20b) 


(20c) 


As  I  remarked  above,  in  the  case  where  additive  noise  is  present  we  would  like  to  regularize  the  expansion 
Eq.(20c)  by  setting  a  threshold  e  >  0  and  using  only  singular  values  and  associated  singular  vectors  for  which 
op  >  e.  The  sensitivity  of  the  inversion  O  to  additive  noise  or  other  errors  in  the  data  vector  is  apparent 
from  Eq.(20c)  which  shows  that  small  changes  in  <  up,  D  >cn  are  magnified  by  the  presence  of  op  in  the 
denominator  of  the  expression  for  O.  Such  a  regularization,  which  is  data  independent,  not  only  increases 
the  robustness  of  the  reconstruction  but  it  also  decreases  the  number  of  singular  functions  that  need  to  be 
stored  so  that  both  computer  storage  requirements  as  well  as  CPU  execution  time  are  decreased.  We  will 
examine  the  regularization  issue  for  a  number  of  test  cases  in  section  6. 


6  Computer  Simulation 

The  inverse  scattering  algorithms  developed  in  the  grant  and  outlined  above  were  tested  in  computer  simu¬ 
lations  in  both  [1]  and  [3].  Here  we  will  only  present  the  simulations  given  in  [3]  since  they  encompass  the 
results  presented  in  the  earlier  paper  and  are  generated  from  the  optimal  form  of  the  reconstruction  algo¬ 
rithm.  In  the  simulations  we  employed  a  simple  non-dispersive  background  model  consisting  of  a  uniform, 
constant  velocity  medium  with  and  without  a  perfectly  reflecting  bottom  layer  in  which  is  embedded  a  small 
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non-dispersive  scattering  potential.  The  scatterer  is  assumed  to  be  embedded  in  a  rectangular  grid  with  x 
denoting  horizontal  position  and  z  vertical  position  with  positive  z  directed  downward.  In  cases  where  we 
employ  a  set  of  frequencies  w/,  /  =  1,2,---  ,  Nu  we  define  a  center  frequency  uic  which  results  in  a  center 
wavelength  Ac  =  1  from  which  we  will  reference  all  of  our  distances.  In  particular,  the  grid  spacing  is  selected 
to  be  one  tenth  of  this  wavelength  (<5x  =  6z  =  1/10).  In  all  simulations  the  transmitter  array  was  located 
on  the  top  row  of  the  rectangular  grid  and  the  receiver  array  along  the  left  hand  side  of  this  grid.  When 
present,  the  reflecting  bottom  was  placed  along  the  bottom  row  of  the  grid. 

The  object  of  the  simulations  is  first  to  examine  the  performance  of  the  algorithm  for  frequency  inde¬ 
pendent  scattering  potentials  when  multiple  distinct  frequencies  are  used  as  discussed  in  section  4  and  to 
compare  the  reconstructions  obtained  in  this  case  with  results  obtained  using  only  a  single  frequency.  A  sec¬ 
ond  goal  is  to  determine  the  computer  storage  gain  obtained  when  using  the  SVD  based  inversion  algorithm 
of  section  5.  In  all  the  simulations  we  assume  point  source  excitation  and  point  receivers  and  do  not  test  the 
effect  of  extended  sources  or  receivers  on  the  performance  of  the  inversion  algorithm.  Also,  as  was  the  case 
in  [1],  the  simulations  test  the  performance  of  the  inversion  algorithm  and  not  the  validity  of  the  DWBA 
which  is  outside  the  scope  of  this  paper.  We  therefore  computed  the  scattered  field  data  (multistatic  data 
matrix)  within  the  DWBA  according  to  Eq.(l)  using  the  above  mentioned  background  consisting  of  either  a 
uniform  constant  velocity  background  or  a  uniform  background  with  a  bottom  reflecting  layer. 


6.1  Test  Cases 


The  simulations  were  performed  on  a  5AC  x  6AC  image  grid  (=  V)  with  one-tenth  center  wavelength  sampling 
as  discussed  above.  The  sources  were  located  on  the  top  row  of  the  support  grid  with  cases  of  Na  —  [3, 5, 7, 9] 
elements  while  the  receivers  we  located  on  the  left  hand  side  of  the  grid  also  with  cases  Np  =  [3,5, 7,9]. 
The  simulations  were  also  performed  at  a  number  of  frequencies  ui  =  [2wc,  wc,  |wc,  Awc]  which  result  in 
wavelengths  of  A  =  [|  Ac,  Ac,  |AC,  2AC]  with  Ac  set  to  a  value  of  unity. 

We  used  two  objects  with  different  characteristics.  First  was  a  composite  object  consisting  of  a  set  of 
5  real  valued  circular  discs  having  diameter  Ac/2  distributed  at  the  corners  and  center  of  a  square  3AC  per 
side.  We  will  call  this  Object  1.  The  second  object  is  a  rectangle  centered  in  V  with  dimensions  3AC  x  2AC. 
Embedded  in  this  square  is  another  rectangle  of  dimensions  1AC  x  Ac/2.  We  will  call  this  Object  2.  Tests 
for  both  objects  were  performed  both  with  and  without  the  reflecting  bottom  and  also  with  and  without 
additive  white  Gaussian  measurement  noise  (AWGN)  although  we  only  present  the  results  for  the  cases  that 
provide  information  we  did  not  see  in  [1]. 

As  was  the  case  in  [1]  the  Green  functions  needed  to  compute  the  spanning  functions  {7rn(r),  n  = 
1,2,  •  ■  ■  ,N}  are  constructed  from  a  Hankel  function  of  the  first  kind  whose  argument  is  a  function  only  of 
the  distance  between  the  field  point  (or  its  mirror  image  about  the  reflecting  bottom)  and  the  source  or 
receiver  location.  In  the  case  of  a  uniform  background  model  the  spanning  functions  are  found  to  be  given 
by 


7Tn(r,W/) 


Gb(Pj  -  r,  w/)G5(afc  -  r,  w/)  if  r  G  V 
0  else-wise 


(21) 


where 

G0{R,uf)  =  ~H0(kUf\R\) 

where  H0  is  the  zero  order  Hankel  function  of  the  first  kind  and  kWj  is  the  background  wavenumber  for  the 
given  frequency  w/ .  In  the  case  of  a  perfect  reflecting  bottom  the  background  Green  function  is  given  by 


Gb(r,r',wf)  =  G0(|r-r'|,w/)  -  G0(|r  -  r'|,  w/) 


where  we  have  used  the  subscript  “b”  to  denote  the  reflecting  bottom  Green  function  and  where  r'  is  the 
mirror  image  of  the  point  r'  about  the  bottom;  i.e.,  r'  =  (x',  —z'  -  2 10)  where  lo  is  the  distance  of  the  bottom 
from  the  top  of  the  grid.  In  the  case  of  a  reflecting  bottom  the  spanning  functions  axe  thus  computed  using 
Eq.(21)  with  the  background  Green  function  replaced  by  G&. 

To  analyze  the  performance  of  the  algorithm  in  the  presence  of  noise,  we  used  additive  white  Gaussian 
noise  (AWGN)  which  was  added  to  the  data  vector  D  such  that  the  noisy  data  Dv  =  D  +  J\f  and  A f  = 
Af(0,  a2 1).  The  noise  variance  was  taken  to  be  2%  of  the  largest  component  of  the  D  vector. 
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6.2  Multiple  Frequencies 

The  use  of  multiple  frequencies  in  the  reconstruction  of  non-dispersive  objects  is  illustrated  in  figures  1-3.  In 
figures  1  and  2,  we  see  the  difference  between  the  reconstructions  for  Object  1  and  Object  2  when  we  use  the 
data  obtained  by  performing  the  scattering  experiment  with  one  frequency  and  then  with  four  frequencies. 
As  we  would  expect,  the  reconstruction  quality  improves  dramatically  as  we  add  the  new  data.  In  particular, 
the  extended  square,  Object  2,  which  is  completely  unresolved  at  a  single  frequency,  is  well  resolved  with 
four.  An  interesting  thing  to  note  is  the  improvement  in  the  ability  of  the  algorithm  to  represent  the  part 
of  the  object  that  is  farther  away  from  the  array.  This  is  critical  as  most  practical  applications  will  not 
have  array  coverage  on  all  sides  of  the  image  grid  and,  as  we  saw  from  part  I,  even  large  arrays  at  a  single 
frequency  will  not  generate  adequate  data  to  accurately  recreate  the  far  edge  of  the  grid.  These  results  are 
consistent  with  other  attempts  at  using  multiple  frequency  data  within  the  DWBA  [13]. 

In  section  4  we  mentioned  that  in  reflection  type  experiments  frequency  diversity  is  more  important 
than  source  receiver  diversity.  The  offset  VSP  type  of  experiment  that  we  are  simulating  here  is  inter¬ 
mediate  between  a  reflection  and  transmission  type  of  experiment  [8];  when  we  use  the  reflecting  bottom, 
the  experiment  is  more  like  a  reflection  experiment  and  without  the  bottom  it  is  more  like  a  transmission 
experiment.  We  should  then  expect  that  frequency  diversity  is  a  more  important  factor  in  the  quality  of  the 
reconstruction  when  the  reflecting  bottom  is  used  and  diversity  in  the  source/receiver  positions  to  be  more 
important  without  the  bottom.  To  illustrate  this,  we  performed  simulations  with  the  reflecting  bottom  and 
with  comparable  amounts  of  data,  one  diverse  in  frequency  and  the  other  diverse  in  receiver  position.  The 
results  for  the  reflection  type  experiment  are  shown  in  figure  3  while  the  results  of  the  transmission  type 
experiment  are  shown  in  figure  4.  In  figure  3(a,b)  we  can  see  the  reconstructions  for  a  7  x  7  array  at  a  single 
frequency  and  for  a  5  x  5  array  at  two  frequencies  which  use  49  and  50  data  points  respectively.  In  figure 
3(c,d)  we  used  a  9  x  9  single  frequency  array  and  a  5  x  5  array  at  three  frequencies,  which  use  81  and  75 
data  points.  After  comparing  these  reconstructions,  we  can  see  that  the  diverse  frequency  arrays  outperform 
the  single  frequency  arrays.  The  single  frequency  arrays  do  a  better  job  of  representing  the  disk  in  the 
upper  left  corner  but  cannot  give  any  indication  of  the  two  discs  along  the  bottom.  In  figure  4  we  do  the 
same  comparison  without  the  reflecting  bottom.  It  is  clear  from  this  figure  that  the  diversity  in  the  position 
of  sources/receivers  produce  a  much  sharper  picture  in  this  experiment  which  again  is  closely  related  to  a 
transmission  experiment. 

6.3  SVD  Inversion  Scheme 

To  examine  the  performance  of  the  SVD  based  formulation  of  the  inversion  scheme,  we  investigated  the 
consequences  of  changing  the  level  of  e.  For  noise-free  data,  as  we  set  the  threshold  higher  and  higher,  the 
number  of  singular  values  we  keep  decreases  and  therefore  the  quality  of  the  reconstruction  will  decrease. 
At  the  same  time  however,  the  time  to  compute  the  reconstruction  as  well  as  the  amount  of  storage  needed 
for  the  singular  functions  decreases,  which  is  the  goal  of  the  new  formulation.  In  fact,  the  storage  space  and 
reconstruction  time  are  both  linear  functions  of  the  number  of  singular  values  with  the  same  rate  of  change, 
so  we  can  think  of  these  qualities  as  equivalent  in  terms  of  their  change.  In  the  case  of  noisy  data,  as  we 
set  the  threshold  lower  and  lower,  we  will  get  a  better  reconstruction  until  a  point  when  the  singular  values 
become  too  small  and  the  reconstruction  is  dominated  by  noise.  It  is  at  this  point  that  we  would  like  to 
compute  the  reconstructions  and  determine  how  much  storage  space  was  saved. 

We  first  performed  simulations  with  noise-free  data,  varying  e  so  that  in  each  test  10%  of  the  singular 
values  were  dropped.  Looking  at  figures  5  and  6,  we  can  see  the  results  of  these  tests.  We  see  that  the 
reconstruction  is  almost  unchanged  even  when  we  remove  50%  of  the  singular  values.  While  it  is  obvious 
that  the  quality  of  the  reconstruction  will  drop  as  we  remove  singular  values,  we  are  also  interested  in  how 
this  compares  to  the  storage  space  needed  for  the  singular  functions.  Looking  at  figure  7  we  can  see  these 
results  for  our  tests.  Here  we  have  plotted  the  percentage  decrease  in  the  time  to  compute  the  reconstruction 
for  some  7V%  of  the  singular  values  defined  as 


<5T/v%  = 


Tiqq%  -  Tn% 


where  Tjoo%  is  the  time  compute  the  reconstruction  using  all  singular  values,  against  the  percent  reconstruc- 
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tion  error,  defined  as 


Er 


lQ(r)-6(r)l2 

|0(r)P 


5x8  z 


]Q(rra)-0(rn)l2 

|0(r„)|2 


(22) 


We  see  that  for  each  of  the  tests,  as  we  remove  10%  of  the  singular  values  for  a  test,  we  also  get  a  10% 
reduction  in  the  time  to  compute  the  reconstruction.  What  is  interesting  is  that  the  percentage  error  is  not 
increasing  at  that  fast  of  a  rate,  especially  in  the  range  of  50  —  100%  of  singular  values.  This  tells  us  that 
when  keeping  between  50  —  100%  of  the  singular  values,  for  a  reasonably  small  decrease  in  reconstruction 
quality,  we  get  a  very  large  performance  increase  in  terms  of  reconstruction  time. 

In  the  noisy  cases,  we  need  to  determine  the  threshold  at  which  we  need  to  start  to  drop  singular  values 
to  obtain  a  reasonable  reconstruction.  This  would  give  us  a  clear  picture  of  how  much  storage  and  time  we 
could  save  while  still  producing  the  optimal  reconstruction.  In  figure  8,  we  can  see  the  noisy  reconstructions 
of  Object  2  for  four  sets  of  frequencies,  two  of  which  we  have  seen  the  in  the  noise-free  cases  of  figure  1.  We 
can  see  here  that  the  objects  were  still  reconstructed  quite  well,  although  we  also  see  that  there  is  a  limit 
to  how  much  data  we  can  add  for  a  given  experiment  as  the  addition  of  the  fourth  frequency  in  figure  8(d) 
does  not  add  much  to  the  reconstruction.  The  main  purpose  of  these  experiments  was  to  determine  how 
many  of  the  singular  values  we  are  able  to  discard,  and  by  looking  at  figure  9,  we  can  see  that  we  were  able 
to  eliminate  5  —  30%  of  the  singular  values,  corresponding  to  an  equivalent  reduction  in  computation  time 
and  storage  space.  We  can  also  see  that  the  percentage  kept  was  well  within  the  region  we  defined  above 
as  having  a  small  impact  on  reconstruction  quality  which  is  why  we  were  able  to  get  quality  reconstructions 
while  reducing  the  storage  needs  by  this  amount. 

It  is  also  important  to  mention  that  the  thresholds  that  were  obtained  were  done  so  with  the  use  of  the 
object  profile  which  allows  us  to  calculate  an  error  value  for  each  amount  of  singular  values  kept.  Without 
knowledge  of  the  object  profile,  the  calculation  of  the  threshold  becomes  a  much  harder  problem  and  is 
beyond  the  scope  of  this  paper.  Attempts  to  work  on  this  problem  are  detailed  in  [14] . 


7  Summary 

In  this  grant  the  inverse  scattering  problem  was  formulated  and  solved  within  the  distorted  wave  Born 
approximation  (DWBA).  The  theory  is  general  and  applies  to  both  dense  and  sparse  transmitter  and  receiver 
arrays  that  are  arbitrarily  distributed  in  space  and  whose  elements  can  be  extended  and  have  arbitrary 
radiation  patterns.  Inversion  algorithms  have  been  derived  and  tested  in  computer  simulations  that  yield  a 
pseudo-inverse  of  the  inverse  scattering  problem  from  either  single  frequency  or  broad  band  scattering  data. 
An  important  feature  of  the  theory  and  inversion  algorithms  is  that  they  incorporate  a  priori  knowledge 
of  both  the  background  in  which  the  unknown  scatterers  are  embedded  as  well  as  estimates  of  the  support 
volume  of  the  scatterers. 

Future  work  in  this  area  will  focus  on  its  extension  to  vector  valued  fields  (e.g.,  electromagnetic  and 
elastic  waves)  and  to  the  use  of  regularization  methods  other  than  the  SVD.  It  is  also  planned  to  apply  the 
algorithm  to  real  data  sets  acquired  in  both  electromagnetic  as  well  as  elastic  wave  scattering  experiments 
currently  being  planned. 
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Figure  1:  Reconstructions  of  Object  1  for  multiple  frequencies  with  reflecting  bottom.  Shown  are  surface  plots 
of  the  absolute  value  of  the  reconstructions  for  (a)  A  =  [Ac]  (b)  A  =  [|AC,  Ac,  |  Ac,  2AC].  Both  reconstructions 
performed  for  Na  =  Np  =  9. 
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Horizontal  Position  (X/10) 


Figure  2:  Reconstructions  of  Object  2  for  multiple  frequencies  with  reflecting  bottom.  Shown  are  surface  plots 
of  the  absolute  value  of  the  reconstructions  for  (a)  A  =  [Ac]  (b)  A  =  [|  Ac,  Ac,  |AC,  2AC] .  Both  reconstructions 
performed  for  Na  =  =  9. 
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Figure  3:  Reconstructions  of  Object  1  comparing  single  and  multiple  frequencies  with  reflecting  bottom. 
Shown  are  the  reconstructions  for  (a)  Na  =  Np  =  7  and  A  =  [Ac]  (b)  Na  —  Np  =  5  and  A  =  [Ac,  |AC]  (c) 
Na  =  Np  =  9  and  A  =  [Ac]  (d)  Na  =  Np  =  5  and  A  -  [|AC,  Ac,  §AC]. 
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Figure  4:  Reconstructions  of  Object  1  comparing  single  and  multiple  frequencies  with  no  reflecting  bottom. 
Shown  are  the  reconstructions  for  (a)  Na  —  Np  =  7  and  A  =  [Ac]  (b)  Na  =  Np  =  5  and  A  =  [Ac,  |AC]  (c) 
Na  =  Np  =  9  and  A  =  [Ac]  (d)  Na  =  Np  -  5  and  A  =  [|AC,  Ac,  §AC]. 


15 


a 


Horizontal  Position  (U10) 


b 


Horizontal  Position  (W10) 


Figure  5:  Surface  plots  of  the  reconstruction  of  Object  1  for  (a)  50%  and  (b)  100%  of  singular  values  kept 
and  noise  free  data.  Reconstructions  performed  with  Na  =  Np  —  7  and  A  =  [|AC,  Ac,  |AC]. 
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Figure  6:  Reconstructions  of  Object  2  for  varying  number  of  singular  values  kept  and  noise  free  data. 
Shown  are  the  reconstructions  for  (a)  25%  (b)  50%  (c)  75%  (d)  100%.  All  reconstructions  performed  for 
Na  =  Np  =  7  and  A  =  [±AC,  ACI  §AC]. 
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Figure  7:  Graph  of  percentage  time  gained  (5T^ r%)  and  percentage  reconstruction  error  (E r)  for  varying 
numbers  of  singular  values  kept.  All  data  was  obtained  for  reconstructions  using  Na  =  Np  =  7  and 
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Figure  8:  Reconstructions  of  Object  1  for  multiple  frequencies  with  noisy  data  and  reflecting  bottom.  Shown 
are  the  reconstructions  for  (a)  A  =  [Ac]  (b)  A  =  [Ac,  |AC]  (c)  A  =  [|AC,  Ac,  | Ac]  (d)  A  =  [^Ac,  Ac,  |AC)  2Ac]. 
All  reconstructions  performed  for  Na  —  Np  =  5. 


Figure  9:  Graph  of  the  percentage  of  singular  values  kept  for  the  reconstructions  of  Object  2  as  seen  in  Fig. 
8  plotted  against  the  number  of  total  singular  values. 
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